# Dependencies ----
library(sf)
library(tidyverse)
library(modelsummary)
library(texreg)
library(gridExtra)

# Load datasets -----
parades <- get(load(file = "data/ICF_dataset")) # ICF dataset
bes_survey <- get(load(file = "data/bes_survey")) # British election Study
beyond_conflict_survey <- get(load(file = "data/beyond_conflict_survey")) # 2023 survey

# Shape files for visualization
ni <- read_sf("data/shapefiles/OSNI_Open_Data_-_Largescale_Boundaries_-_NI_Outline//") %>%
  st_make_valid()
constituencies <- st_read("data/shapefiles/OSNI_Open_Data_-_50K_Boundaries_-_Parliamentary_Constituencies/") %>%
  mutate(constituency = tolower(PC_NAME)) %>%
  mutate(constituency = case_when(
    constituency == "down south" ~ "south down",
    constituency == "down north" ~ "north down",
    constituency == "newry and armagh" ~ "newry & armagh",
    constituency == "fermanagh and south tyrone" ~ "fermanagh & south tyrone",
    TRUE ~ constituency)) %>%
  st_make_valid()



# Replication code ----
## Figure 1 ----
df_processed <- parades %>%
  count(year, protestant) %>%
  mutate(parade_type = ifelse(is.na(protestant), "Total parades", 
                              ifelse(protestant == F, "Catholic parades", "Protestant parades"))) %>%
  bind_rows(
    parades %>%
      count(year, sensitive) %>%
      filter(sensitive == TRUE) %>%
      mutate(parade_type = "Sensitive parades") 
  ) %>%
  bind_rows(
    parades %>%
      count(year) %>%
      mutate(parade_type = "Total parades")
  )

df_processed$parade_type <- factor(df_processed$parade_type, 
                                   levels = c("Total parades", "Sensitive parades", "Catholic parades", "Protestant parades"))

figure1 <- ggplot(df_processed) + 
  geom_line(aes(x = year, y = n, group = parade_type, linetype = parade_type, color = parade_type), linewidth = 1) +
  scale_linetype_manual(values = c("Total parades" = 1, "Sensitive parades" = 2, "Catholic parades" = 3, "Protestant parades" = 6)) +  
  scale_color_grey() +  
  scale_x_continuous(name = "", breaks = seq(2001, 2022, by = 2)) +
  scale_y_continuous(name = "Number of parades", limits = c(0, 3400), breaks = seq(0, 3200, by = 400)) +
  theme_classic() +
  theme(legend.title = element_blank(),
        legend.position = "bottom",
        text = element_text(size = 15))


## Figure 2 ----
figure2 <- parades %>%
  st_as_sf() %>%
  st_transform(crs = st_crs(ni)) %>%
  filter(st_intersects(ni, ., sparse = FALSE)[1,]) %>%
  ggplot() +
  geom_sf(color = "black", alpha = .05) +
  geom_sf(data = ni, fill = NA) +
  theme_void() +
  theme(legend.position = "bottom")
figure2


## Figure 3 ----
# Combine parades and constituencies
parades_constit <- constituencies %>%
  st_make_valid() %>%
  st_join(parades, join = st_intersects) %>%
  st_drop_geometry() %>%
  group_by(constituency, year) %>%
  summarise(
    num_parades = n(),
    sensitive = mean(sensitive, na.rm = TRUE),
    counter_protest = mean(counter_protest, na.rm = TRUE),
    conditions_imposed = mean(conditions_imposed, na.rm = TRUE),
    protestant = mean(protestant, na.rm = TRUE)
  )

sensitive_parades_constit <- constituencies %>%
  st_make_valid() %>%
  st_join(parades, join = st_intersects) %>%  
  st_drop_geometry() %>%
  filter(sensitive == TRUE) %>%
  group_by(constituency, year) %>%
  summarise(
    num_sensitive_parades = n(),
    counter_protest = mean(counter_protest, na.rm = TRUE),
    conditions_imposed = mean(conditions_imposed, na.rm = TRUE),
    protestant = mean(protestant, na.rm = TRUE)
  ) %>%
  left_join(parades_constit %>% select(constituency, year, num_parades), by = c("constituency", "year")) %>%
  mutate(prop_sensitive = ifelse(num_parades > 0, num_sensitive_parades / num_parades, 0)) %>%
  ungroup()


figure3 <- constituencies %>% 
  merge(data.frame(year = unique(bes_survey$year)), by = NULL) %>%
  left_join(sensitive_parades_constit %>% filter(year %in% unique(bes_survey$year))) %>%
  mutate(num_sensitive_parades = replace_na(num_sensitive_parades, 0)) %>%
  ggplot() +
  geom_sf(aes(fill = num_sensitive_parades), color = "black") +
  facet_wrap(~year, nrow = 3) +
  theme_void() + 
  scale_fill_gradient(
    low = "white", high = "black", na.value = "grey50",
    name = "Contentious parades\nper constituency",
    limits = c(0, 100),
    breaks = seq(0, 100, 25), labels = c("0", "25", "50", "75", ">100")
  ) +
  theme(legend.position = 'bottom') +
  guides(fill = guide_colourbar(title.position="top", title.hjust = 0.5))



## Table 2 ----
m1 <- lm(mixed_schooling ~  num_sensitive_parades + num_parades + age + female + education + religion + year_f, bes_survey)
m2 <- lm(marraige_relative ~ num_sensitive_parades + num_parades + age + female + education + religion + year_f, bes_survey)

texreg::screenreg(list(m1, m2), 
                  custom.coef.map = list(
                    "num_sensitive_parades"  = "Num. contentious parades (in 100s)",
                    "prop_sensitive"  = "Prop. contentious parades"
                  ), 
                  stars = c(0.001, 0.01, 0.05, 0.1))


## Figure 4 ----
m1 <- lm(mixed_schooling ~  num_sensitive_parades + num_parades + age + female + education + religion, bes_survey %>% filter(year==2003))
m2 <- lm(mixed_schooling ~  num_sensitive_parades + num_parades + age + female + education + religion, bes_survey %>% filter(year==2010))
m3 <- lm(mixed_schooling ~  num_sensitive_parades + num_parades + age + female + education + religion, bes_survey %>% filter(year==2015))
m4 <- lm(mixed_schooling ~  num_sensitive_parades + num_parades + age + female + education + religion, bes_survey %>% filter(year==2017))
m5 <- lm(mixed_schooling ~  num_sensitive_parades + num_parades + age + female + education + religion, bes_survey %>% filter(year==2019))

df1 <- modelsummary::modelplot(list("2003" = m1,
                                    "2010" = m2,
                                    "2015" = m3,
                                    "2017" = m4,
                                    "2019" = m5),
                               coef_map = c(num_sensitive_parades = "Num. contentious parades (in 100s)"), draw = F)


m6 <- lm(marraige_relative ~  num_sensitive_parades + num_parades + age + female + education + religion, bes_survey %>% filter(year==2003))
m7 <- lm(marraige_relative ~  num_sensitive_parades + num_parades + age + female + education + religion, bes_survey %>% filter(year==2010))
m8 <- lm(marraige_relative ~  num_sensitive_parades + num_parades + age + female + education + religion, bes_survey %>% filter(year==2015))
m9 <- lm(marraige_relative ~  num_sensitive_parades + num_parades + age + female + education + religion, bes_survey %>% filter(year==2017))
m10 <- lm(marraige_relative ~  num_sensitive_parades + num_parades + age + female + education + religion, bes_survey %>% filter(year==2019))

df2 <- modelsummary::modelplot(list("2003" = m6,
                                    "2010" = m7,
                                    "2015" = m8,
                                    "2017" = m9,
                                    "2019" = m10),
                               coef_map = c(num_sensitive_parades = "Num. contentious parades (in 100s)"), draw = F)


figure4 <- rbind (df1 %>% mutate(DV = "Mixed schooling"),
             df2 %>% mutate(DV = "Mixed marriage")) %>%
  ggplot(aes(model, estimate)) +
  geom_point(size = .1) +
  facet_wrap(~DV, nrow = 2) +
  geom_linerange(aes(ymin = conf.low, ymax = conf.high), linewidth = .2) +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_x_discrete(name = "") +
  scale_y_continuous(name = "") +
  theme_classic() +
  theme(legend.position = "none")



## Figure 5 ----

p1 <- ggplot(beyond_conflict_survey) +
  geom_jitter(aes(survey_date, y = 0, group = Case, color = surveyed_after), height = 0.3, size = .5) +
  geom_vline(xintercept = dmy("12/07/2022")) + 
  scale_color_manual(values = c("black", "grey")) +
  scale_y_continuous(name = "", limits = c(-.5, .5)) +
  scale_x_date(labels = scales::date_format("%d %b"), date_breaks = "2 weeks") +
  labs(x = "") +
  theme_classic() +
  theme(legend.position = "non", 
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

p2 <- ggplot(beyond_conflict_survey, aes(x = factor(surveyed_after, levels = c(F, T)), y = mixed_living-1, color = surveyed_after)) +
  stat_summary(fun.data = "mean_se") +
  scale_color_manual(values = c("black", "grey")) +
  #coord_flip() +
  scale_x_discrete(name = "", labels = c(
    "Interviewed before\nthe Twelfth",
    "Interviewed after\nthe Twelfth")) +
  scale_y_continuous(name = "Mixed neighbourhood", labels = scales::percent_format(accuracy = 1),
                     breaks = seq(from = 0.45, to =  .8, by= .015)) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text.y = element_text(hjust = .5),
        panel.grid.minor = element_line(linewidth = 0.25, linetype = 'solid',
                                        colour = "grey90"))

gridExtra::grid.arrange(p1,p2, nrow = 2)


## Figure 6 ----
result_list = list()
for (i in c(1:8)) {
  df <- beyond_conflict_survey %>% 
    filter((survey_date > (dmy("12/07/2022")-weeks(4))) & (survey_date < dmy("12/07/2022")+ weeks(i)))
  
  m1 <- lm(mixed_living ~ surveyed_after + age + factor(gender) + income + education, df)
  m2 <- lm(mixed_living ~ surveyed_after + age + factor(gender) + income + education, df %>% filter(religion=="Catholic"))
  m3 <- lm(mixed_living ~ surveyed_after + age + factor(gender) + income + education, df %>% filter(religion=="Protestant"))
  
  
  results1 <- data.frame(
    estimate = coef(m1),
    std_error = summary(m1)$coefficients[, "Std. Error"])
  results1$variable = row.names(results1)
  results1$model = "All respondents"
  
  results2 <- data.frame(
    estimate = coef(m2),
    std_error = summary(m2)$coefficients[, "Std. Error"])
  results2$variable = row.names(results2)
  results2$model = "Catholic respondents"
  
  results3 <- data.frame(
    estimate = coef(m3),
    std_error = summary(m3)$coefficients[, "Std. Error"])
  results3$variable = row.names(results3)
  results3$model = "Protestant respondents"
  
  
  results <- rbind(results1, results2, results3) %>% cbind(temporal_window = rep(i, nrow(results1)))
  
  result_list[[i]] = results
  
}

figure6 <- do.call(rbind, result_list) %>%
  filter(variable %in% c("surveyed_afterTRUE")) %>%
  mutate(variable = ifelse(variable=="(Intercept)", "Intercept", "Surveyed after the Twelfth")) %>%
  ggplot(aes(temporal_window, estimate, color = model)) +
  geom_point() +
  geom_linerange(aes(ymin=estimate+std_error*2, ymax=estimate-std_error*2)) +
  facet_wrap(~model, scales = "free_y") +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_y_continuous(limits = c(-.5, .5)) +
  scale_x_continuous(name = "Number of weeks after 12th July",
                     breaks = c(1:9)) +
  scale_color_grey(name = " ") +
  theme_classic() +
  theme(legend.position = "none", strip.background = element_blank(),
        strip.text.x = element_text(face = "bold"))

figure6


